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tional integrals each giving the contribution from quantum states with definite symmetry proper- 
ties. We propose a strategy to compute each of them, normalized to the one with vacuum quantum 
numbers, by a Monte Carlo procedure whose cost increases power-like with the time extent of the 
lattice. This is achieved thanks to a multi-level integration scheme, inspired by the transfer matrix 
formalism, which exploits the symmetry and the locality in time of the underlying statistical sys- 
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Carlo. We briefly illustrate the approach in the simple case of the one-dimensional harmonic os- 
cillator and discuss in some detail its extension to the four-dimensional Yang Mills theories. We 
report on our recent new results in the SU(3) Yang-Mills theory on the relative contribution to the 
partition function of the parity-odd states. 
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1. Introduction 

1.1 The problem 

Lattice field theories can be studied numerically by Monte Carlo simulations. They allow to address 
non-perturbative problems from first principles, and for most of the theories the lattice provides the 
only known non-perturbative definition. The mass of the lowest states in a given channel can, 
for instance, be extracted from the Euclidean time-dependence of a suitable two-point correlation 
function. Very often, however, the statistical error of the Monte Carlo estimate grows exponen- 
tially with time, and in practice it is not possible to find a window where statistical and systematic 
errors are both under control. The problem is easily explained by looking at the states contributing 
asymptotically in time to the two-point function and to its variance. Whenever the energy of the 
asymptotic state in the variance is smaller than twice that in the two-point function, the noise to 
signal ratio is going to grow exponentially in time M, g]. The standard Monte Carlo approach fails 
basically because for any given field configuration all asymptotic states of the theory are allowed 
to propagate in the time direction, regardless of the quantum numbers of the source fields. Their 
contributions disappear in the Monte Carlo average for the two-point function but sum up in the 
noise. As shown in the following subsection, the issue is already there for a simple system such as 
the harmonic oscillator. We use the latter to introduce the basic ideas of the method that was pro- 
posed in [Q, H], a "symmetry-constrained" Monte Carlo, and to show how it avoids the exponential 
increase of the signal-to-noise ratio. 

1.2 The case of the harmonic oscillator 

We consider the one-dimensional harmonic oscillator on the lattice. We recall here a few basic 
equations. More details can be found in [^] to which we refer for any unexplained notation. The 
system is described by the Hamiltonian 

H = — + V(x) with V(x) = -m(D 2 x 2 . (1.1) 
2m w w 2 

This operator is invariant under parity transformations, therefore its eigenstates can be classified 
according to a parity quantum number (+ or — ). We label the corresponding energy levels as 
and Sj respectively. The transfer operator between two consecutive time slices is defined as 

1 # = g -|v(*) e -«i e -Sv(D > (12) 
an its matrix elements in the coordinate basis 

/ m \ 1/2 

(x n+l \#\x n ) = J T n+hn (1.3) 

can be computed explicitly 

T n+ i,„ = e- aL »^» , (1.4) 



with 

L n+ , n s S?( Xn+uXn ) = % ( 2 + ^ + ^ • (1.5) 

2 \ a 2 2 
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The statistical variance associated to the two-point correlation function {xix k ) (interpolating 
parity odd states) is 

a 2 = (xf4) - (x lXk ) 2 , (1.6) 

and, at asymptotically large time separations, the signal-to-noise ratio can be easily computed in 
the underlying quantum field theory 

( Xl x k ) _ |(4"|i|4 + )| 2 (# -_^ ) | / _^| 

a |<4Wo + >l 

The ratio decreases exponentially in time, as announced. As we will describe in the following the 
problem here can be solved by introducing the "sign" and the "module" fields. 
We define the complete set of parity eigenstates 

\x,±) = ^=(\x)±\-x}), &\x,±) = ±\x,±). (1.8) 

The invariance of the Hamiltonian under parity implies 

/2m \ l / 2 

(s',x n+l \^\x n ,s) = i^—j T° +hn 8 s , s , (1.9) 



with 



A n+l,n = \{^(-X„ + l,X n )±^(x n+l ,X n )} , (1.11) 



and the functional integral can be written as 



s=± J n=0 



where N is the extent of the lattice. We further define 

-aL+ , „ _ -aA" 



T n \i,n = = e~ aA " +l " cosh{«A„- +1 „} , (1.13) 

and cast the functional integrals in the form 



t.N-1 r N-l N-l 

1+ = / n dx ne- S+ , Z = U dx n e- s+ \{ tanh{a A m+1 J , (1.14) 

J n=0 J n=0 m=0 

7+ = z,^- 1 7" + 



where S + = aE„=o %+i n- ^ ne P atn i nte g ra l is thus rewritten as a sum of two functional integrals 
giving the contribution from parity even and odd states respectively. Each integrand is a product 
of transfer matrix elements between quantum states with definite parity. The two-point correlation 
function reads (k < I) 

l r N ~ l + r /_1 
(xix k ) = - Yl dx n e- s+ \xi Y\twh{aA m+lm }x k 

^ J n=0 k m=k 

N-l k-1 

+ Y[ tanh{«A m+1 m }x[X k Y[ tanh{aA m+1 m }| . (1.15) 



m=l m=0 
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Figure 1: Left: two-point correlation function versus the time distance t/a of the sources. Right: errors 
on the correlation function as obtained with the multi-level algorithm (SCMC) and with the standard path 
integral Monte Carlo (PIMC). 



Each term is now the expectation value of a factorized observable in a system described by the 
action S + . The expression reminds of the factorized form used for the correlator of Polyakov loops 
in the pure gauge theory in As done there a multi-level integration scheme can be introduced 
for the system and the observable discussed here. The key ingredients are sub-lattice averages, i.e. 
averages computed by numerically integrating over the degrees of freedom in a thick time-slice of 
the lattice with the variables at the boundaries kept fixed, and the recursive relations, which allow 
to obtain averages on large thick time-slices as the product of those on smaller ones integrated over 
their boundaries configurations. Both properties are due to the locality of the action. We do not 
repeat here the details concerning the construction of the algorithm, they can be found in [^] but 
rather report on the main results. 

In the left plot of Figure [j] we show the two-point correlation function computed on a lattice 
with N = 64 points, with statistical errors being smaller than symbols. The error (SCMC) is shown 
on the right plot of the same Figure. The signal-to-noise ratio is depleted, as expected, (only) 
inversely proportional to the time distance of the sources. For comparison in the same plot it is also 
shown the statistical error obtained with a standard Monte Carlo procedure (PIMC) which needed 
roughly the same CPU time. It is clear that with our strategy the statistical error is exponentially 
reduced, and at large time distances it is lowered by many orders of magnitude. The effective 
energy-split a(b(t) extracted from the correlator is shown in the left plot of Figure |2j, and it is in 
perfect agreement with the theoretical expectation On the right plot of the same figure it is 
shown the effective estimate R(t) of the square of the matrix element {Sq\x\Sq) computed as 



' T 



R(t) = - 7 F7f 7} 7TTT , d-16) 

2cosh [co (j - a\l - k\) \ 

which also agrees very well with the analytical result. 

The approach described here is of inspiration for systems with a larger number of degrees of 
freedom. In that case, however, one cannot make sense of the "sign" and "module" fields. New 
concepts have to be introduced as we illustrate in the following for the SU (3) Yang-Mills theory 
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Figure 2: Left: effective energy difference extracted from the two-point correlator at any time separation 
t/a. Right: the ratio R as defined in the text. 



2. Extension to Yang-Mills theories on the lattice 

We again divide the parity even sector of the theory from the parity odd, naming the energy levels 
and EJ respectively. The statistical variance of the estimate of a two-point correlation func- 
tion (O(xq)O(O)) of a parity-odd interpolating operator O, computed by the standard Monte Carlo 
procedure, is defined as 

a 2 = (<9 2 (x )O 2 (0)> - (O(x )O(0)} 2 . (2.1) 
At asymptotically-large time separations the signal-to-noise ratio takes the form 



<0(*>)0(o)> _ \(E^\d\o}\\_ E 



l(0|O 2 |0)| 



e ~' X0 + ' 



(2.2) 



i. e. the signal is again depleted exponentially in time. 

For the one-dimensional harmonic oscillator the quantity tanh{aA m+1 m } represented the ratio 
Z~/Z + for a system of one time-slice with fixed boundary configurations. Due to the regularity of 
the spectrum this ratio is of 0(1). The same cannot be expected for systems with many degrees 
of freedom, as the four-dimensional Yang-Mills theory, the ratio will rather be of 0(e~( L / a ' 3 ), with 
L the spatial extent of the lattice. However, if one considers systems d time-slices large, with 
d ~ l/T c and T c the critical temperature, the same ratio is now expected to be of 0(e~ E ' d ) for each 
boundary configuration. These are the quantities we want to directly access and use to rewrite our 
observables. To this end we first need to briefly recall the formalism of the transfer matrix, we refer 
to [Q] for a more thorough discussion. 



2.1 Transfer matrix 

We adopt Wilson's regularization of gauge theories 0. The corresponding transfer matrix has been 
explicitly constructed in |9[ [H], [IT]]. The functional integral with periodic boundary conditions 
in time can be written as 

r-i 



n D 3^o]T 

-vo=0 



V Xo +uV Xo 



(2.3) 
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where the transfer matrix elements among states \V Xo ) in the coordinate basis are defined as 



T 



V, 



x +\,Vx 



with 



















L 


^o+l 5 ^vo 


= K 


V Xo +uV Xo 











(2.4) 



(2.5) 



and \V£) is the result of a gauge transformation Q. on the state \V XQ ). The kinetic and the potential 
contributions to the Lagrangian are given by 



K 



PL l-^ReTr{\4(x +l,x)y;(x ,x)} 



and 



W 



x,k 
P 

2 



LL l-^ReTr{y w (^o,x)} 



x kj 



(2.6) 



(2.7) 



respectively, where Vu is the spatial plaquette computed with the links Vjt(x). By exploiting the 
invariance of the Haar integration measure under left and right multiplication it is easy to show that 
the transfer matrix is gauge invariant. For a thick time-slice, i.e. the ensemble of points in the sub- 
lattice with time coordinates in a given interval [jco^o] and bounded by the equal-time hyper-planes 
at times xq and y^, the transfer matrix elements can be introduced by the formula 



V yo ,V Xo 



3*0-1 



3*0-1 



n o 3 [y W0 ] n t 

wo=*o+l 



zo=*o 



^+1,^0 



The parity transformation acts on the states in the coordinate basis as 

m = |V p ) , |V> = P G \V) , Vf(x) = vi(-x-%) , 



(2.8) 



(2.9) 



where P G is the projector on gauge invariant states. Again, we can then define a complete set of 
parity eigenstates 

1 r. 



|V,±) = — yv)±|v^ 

and their transfer matrix elements are given by 



p\v,±) = ±\y,±) 



t|v^, s > 



(^,V A - +i 



28 s , s r 



■vo 



Vro+l'^o 



(2.10) 

(2.11) 
(2.12) 



For a thick time-slice the matrix elements between parity states can be introduced by exploiting the 
same composition rule as in Eq. (2J5) with T replaced by T s . In addition, the relations 



D 3 [*yT J 



V yo ,V zo 



D 3 [V Z0 ]V V yo ,V a 



o. 



Vy ,V n 



(2.13) 
(2.14) 
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hold. In particular they imply that 

tTv v 1 ~7~ \ D ^ U ^ be Tin tj P ( • 5) 



an useful expression for the practical implementation of the multi-level algorithm described in 
the following. The subscript "sub" indicates that the integral is performed over the dynamical 
field variables in the thick time-slice [jco,;yo] with the spatial components £4(x) of the boundary 
fields fixed to Vk(xo,x) and Vk(yo,x) respectively. Finally, by replacing T^+i,!^,] in Eq. ( jl3"l ) by 



Li^fVio+i, V Xo ] and repeatedly applying Eq. ( 2.13 ), it is possible to rewrite the path integral as a 
sum of functional integrals 



z=£z% z°= /nW^T* 

s=± J x =Q 



(2.16) 



each giving the contribution from gauge-invariant parity-even and -odd states respectively. 

The insertion of T v [V yo , V Xo ] in the path integral plays the role of a projector, as on each config- 
uration it allows the propagation in the time direction of states with parity s only. Indeed the parity 
transformation of one of the boundary fields in T[Vy ,V Xo ] flips the sign of all contributions that it 
receives from the parity-odd states while leaving invariant the rest. The very same applies to the 



path integral in Eq. (2.3) if the periodic boundary conditions are replaced by ^-periodic boundary 



conditions, i.e. Vj = Vq. All contributions from the parity odd states are then multiplied by a 
minus sign. 

2.2 The hierarchical integration scheme 

To determine the parity projector between two boundary fields of a thick time-slice, the basic 
building block to be computed is the ratio of transfer matrix elements 



Wx 0+ d,V X0 ] 



R[V X0+d ,V X0 ] = ^ + " f . (2.17) 



As mentioned above, for d of 0( l/T c ), the ratio R is expected to be of 0( 1). However the integrands 



in the numerator and in the denominator on the r.h.s of Eq. ( J2. 17| ) are, in general, very different 
and the main contributions to their integrals come from different regions of the phase space. The 
most straightforward way for computing R is to define a set of n systems with partition functions 
3?\ ... 3f n designed in such a way that the relevant phase spaces of successive integrals overlap and 
that 3f\ = T \V XQ+ d, VJq] and ^ n = TfV^+^V^,]. The ratio R can then be calculated as 

gp. SP* gp t gp , 
R =^F x ^F x --- x ^p— x ^)^-> ( 2 - 18 ) 

with each ratio on the r.h.s. being computable in a single Monte Carlo simulation by averaging the 
proper reweighting factor. 

For the case at hand 3?\ and % x are the partition functions of two systems differing only for 
the boundary conditions in time. In both cases Dirichlet boundary conditions are imposed but the 
boundary configurations at time xq + d differ by a parity transformation. Instead of relating the 



7 



Symmetries and exponential error reduction in YM theories on the lattice 



boundary configurations in such a way, we change the action of one of the two systems on the 
last time-slice, by introducing a new temporal link connecting the point x,xo +d — 1 on the last 
dynamical time-slice to its parity transformed — x,xo + d on the boundary. We call the associated 
plaquette "parity twisted" space-time plaquette and (parity twisted kinetic term) the sum of 
such plaquettes (see Figure ||). To interpolate between the two systems we slowly switch off the 



k 




Figure 3: Two dimensional representation of the plaquettes in the kinetic terms K (r — 1 /2) and K p (r = 
— 1 /2) on the time-slice yo — 1 . The telescopic algorithm described in the text bridges between the two 
systems in l) steps. 



coupling j8 in K by decreasing it in steps of size e = 1 /L 3 while increasing the coupling in 
by the same amount, we distinguish the interpolating actions by a parameter — 1/2 < r < 1/2. In 
this way we move in L 3 steps from one system to the other. This means we need to perform a 
chain of L 3 Monte Carlo simulations within a Monte Carlo simulation and we therefore have an 
algorithm, which scales with the second power of the volume L 3 . This is the case also for other 
known methods for computing ratios of partition functions p 



13, 14] 



Once the projectors have been computed, the ratio of partition functions Z s /Z can be calculated 
by implementing the hierarchical two-level integration formula 



1 



s 

m,d 



(2.19) 



Where P m,d 



yo,x 



is defined as 



m.d 



yo,xo 



i=0 T [ U x +(i+\)-d,U X0+ i.d} 



(2.20) 



with m > 1 and yo = x o + m ■ d. The procedure can, of course, be generalized to a multi-level 



algorithm. For a three-level one, for instance, each ratio on the r.h.s of Eq. ( |2.20[ ) can be computed 
by a two-level scheme. For each configuration of the boundary fields, the magnitude of the product 
in our observable P^ lrf [r,0] is proportional to T , and the statistical fluctuations are reduced 
to this level. This has to be compared to the standard case in which each configuration gives 
a contribution to the signal which decreases exponentially in time, whereas it contributes 0(1) 
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to the noise (variance) at any time distance. To achieve an analogous exponential gain in the 
computation of the correlation functions, the projectors T s have to be inserted in the proper way 
among the interpolating operators. As a technical aside we remark that the computation of R 
requires a thermalization procedure for each value of r. We do not expect the latter to be particularly 
problematic since, as mentioned above, expectation values for consecutive values of r refer to path 
integrals with the relevant phase spaces which overlap. The ratio R is computed by simulating 
systems corresponding to consecutive values of r one after the other, and by starting from the one 
used to extract the boundary fields (r = 1 /2). 

2.3 Results 

In the four-dimensional SU(3) Yang-Mills theory we have simulated lattices with an inverse gauge 
coupling of j3 = 6/gg = 5.7 which corresponds to a value of the reference scale ro of about 2.93a 
JB] , [T6| ]. The number of lattice points in each spatial direction has been set to L = 6,8 and 10 
corresponding to a linear size of 1.0, 1.4 and 1.7 fm respectively. For each spatial volume we have 
considered several time extents T, the full list is reported in Table [I] together with the number of 
configurations generated, the details of the multi-level simulation algorithm used for each run and 
the results for Z~/Z and for the effective mass M~ of the first parity-odd glueball state extracted 
from that ratio: 

M- = -iln(^(r)V (2.21) 



The natural logarithm of for the interpolating systems in Eq. (|2.18|) is shown as a function 



of r in the left panel of Fig. |] for a typical configuration of the run B3. As expected, its value is of 
0(1) for each value of r. Its almost perfect asymmetry under r — ► — r, however, makes the product 
of all the L 3 results a quantity of 0(1). This impressive cancellation, which is at work for T > 3 on 
all volumes considered, can be better appreciated in the right panel of the same Figure, where the 
sum of the function in the interval [—r,r] is plotted for a subset of values of r. It is the deviation 
from the exact asymmetry which flips in sign under a parity transformation of one of the boundary 
fields, and forms the signal we are interested in. 

The Monte Carlo history of P^ r/ , 2 [r,0] is shown in Figure || for the lattice A5. The central 
dashed line corresponds to the average value, while the other two delimit the one standard devia- 
tion region. As expected the Monte Carlo history is very regular and each configuration gives an 
estimate of the observable which is of the right size. Fluctuations are five times the average value 
at most. We have observed similar Monte Carlo histories also for the other runs. 

Finally we show the results for Z~/Z and aM~ in Figs. ^ and ^ respectively. We have been 
able to follow the exponential decay in the ratio Z~/Z over almost 7 orders of magnitude. The 
data at large values of T/a can be used to estimate the multiplicity of the first parity odd state, a 
quantity which is not accessible within the other approaches. To this end the precision however has 
to be increased, as for now we assume the multiplicity to be one, which justifies the definition of 
the effective mass M~ in Eq. 2.21 . Figure [7] shows that the algorithm works as expected as the error 



on the effective mass could be kept constant to the level of a few percent up to a separation of about 
3.5 fm. It also shows that finite size effects are rather large for lattices of linear size around 1 fm 
(L/a = 6) but they become negligible within the present accuracy once a size of 1 .4 fm (L/a = 8) is 
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Lattice 


L 


T 


A^conf 


Mev 


d 


Z" 

z 


aM~ 


Ai 


6 


4 


50 


2 


4 


0.409(8) 


0.223(5) 


A 2 




5 


50 


2 


5 


0.177(13) 


0.346(14) 


A 3 




6 


50 


2 


6 


0.069(7) 


0.446(17) 


A 4 




8 


175 


2 


4 


1.47(28)-10- 2 


0.528(24) 


A 5 




10 


50 


2 


5 


2.2(5)-10- 3 


0.611(20) 


A 6 




12 


90 


2 


6 


6.6(17)-10- 4 


0.610(21) 


A 7 




16 


48 


2 


8 


2.8(8)-10~ 5 


0.655(18) 


A 8 




20 


48 


3 


{5,10} 


1.5(5)-10 6 


0.670(15) 


Bi 


8 


4 


20 


2 


4 


0.426(8) 


0.213(5) 


B 2 




5 


25 


2 


5 


0.061(6) 


0.558(21) 


B 3 




6 


75 


2 


3 


1.65(26)-10- 2 


0.685(27) 


B 4 




8 


48 


2 


4 


1.37(26)-10- 3 


0.824(24) 


B 5 




12 


48 


3 


{3,6} 


3.6(18)-10 6 


1.045(41) 


B 6 




16 


36 


3 


{4,8} 


5.2(19)-10~ 8 


1.048(23) 


Ci 


10 


4 


20 


2 


4 


0.455(12) 


0.197(6) 


c 2 




5 


24 


2 


5 


0.060(3) 


0.561(11) 


c 3 




6 


50 


2 


3 


1.6(4)- 10- 2 


0.687(39) 


c 4 




8 


48 


2 


4 


5.2(16)-10- 4 


0.944(39) 


c 5 




12 


24 


3 


{3,6} 


3.3(17)10 6 


1.052(43) 



Table 1: Simulation parameters and results. N con f is the number of configurations of the uppermost level, 
Mev is the number of levels and d is the thickness of the thick time-slice used for the various levels. The 
effective mass M~ is given by — T~ l ln(Z~/Z). 

reached. We therefore quote rom G - = 3.07(7) from L/a = 8, T/a = 16 as a preliminary result for 
the mass of the lightest J PC = 0~ + glueball at a lattice resolution of 0.17 fm with Wilson's gauge 
action. Given the quite large value of the lattice spacing, cutoff effects may affect this number 
significantly. 

3. Conclusions and outlook 

For most of the two-point functions computed on the lattice the noise to signal ratio grows expo- 
nentially with the time separation of the source and the sink. This disease can be cured by imposing 
the propagation of states with the desired quantum numbers only on each (gauge) configuration. 
The algorithm proposed here solves the problem by making use of the symmetry properties of 
the underlying quantum theory. We have numerically tested the approach in the four-dimensional 
SU(3) Yang-Mills theory, by computing the mass of the lightest parity-odd glueball. For a given 
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Figure 4: Left: the natural logarithm of is shown as a function of r (statistical errors are smaller than 
symbols) for a typical configuration of the run B3. Right: the sum of the points in the interval [— r,r] is 
plotted as a func 
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Figure 5: Monte Carlo history of the quantity P 2 5 [10, 0] for the run A5 



precision on the latter the algorithm scales as a power of T, the total time extent of the lattice, and 
we have therefore been able to follow an exponential decay over 7 orders of magnitude and up to 
separations of 3.5 fm. That allows to isolate the contribution of a single state with unprecedented 
confidence. We have also studied finite size effect and collected strong indications that, for the 
effective mass considered here and within our statistical errors, those are negligible for L > 1 .4 fm. 
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Figure 6: The quantity Z~/Z vs T/a. 
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Figure 7: The quantity aM vs T/a. 
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The multiplicity of the state can also be computed using the approach described, and in the near 
future we plan to increase the accuracy on its determination, which can be obtained only with lim- 
ited precision by using the data produced so far. The reduction of systematic uncertainties related 
to lattice artifacts remains an expensive task as the algorithm scales roughly as (L/a) 6 . 

The inclusion of other symmetries is straightforward. We have already implemented charge- 
conjugation and tested it in small volumes, observing basically the same efficiency of the integra- 
tion scheme as for the parity discussed here. Different symmetry transformations can be actually 
considered simultaneously and we plan to include cubic rotations and translations. The mass of the 
lightest state in any sector specified by the quantum numbers J PC could then be computed without 
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suffering from the exponential problem JT7| ] . 

A way to generalize the ideas reported here to systems including fermion degrees of freedom 
is, at present, not known. Among other advantages, such an extension would allow to compute 
the ratio between the partition functions at different baryon quantum numbers avoiding the sign 
problem, which affects the simulations at finite density. 

References 

[1] G. Parisi, Phys. Rept. 103 (1984)203. 

[2] G.P Lepage, TASI 89 Summer School, Boulder, CO, Jun 4-30, 1989. 
[3] M. Delia Morte and L. Giusti, Comput. Phys. Commun. 180 (2009) 813. 

[4] M. Delia Morte and L. Giusti, Comput. Phys. Commun. 180 (2009) 819, arXiv:0806.2601 [hep-lat]. 
[5] M. Luscher and P. Weisz, JHEP 09 (2001) 010, hep-lat/0108014. 
[6] M. Creutz and B. Freedman, Ann. Phys. 132 (1981) 427. 
[7] K. G. Wilson, Phys. Rev. D 10, 2445 (1974). 

[8] K.G. Wilson, in "New developments in quantum field theory and statistical mechanics", Cargese 
1976, Eds. M. Levy and P. Mitter, Plenum (NY 1977). 

[9] M. Luscher, Commun. Math. Phys. 54 (1977) 283. 

[10] M. Creutz, Phys. Rev. D15 (1977) 1128. 

[11] K. Osterwalder and E. Seiler, Ann. Phys. 1 10 (1978) 440. 

[12] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195. 

[13] C. Hoelbling, C. Rebbi and V.A. Rubakov, Phys. Rev. D63 (2001) 034506, hep-lat/0003010. 

[14] P. de Forcrand, M. D'Elia and M. Pepe, Phys. Rev. Lett. 86 (2001) 1438, hep-lat/0007034. 

[15] ALPHA Coll., M. Guagnelli, R. Sommer and H. Wittig, Nucl. Phys. B535 (1998) 389, 
hep-lat/9806005. 

[16] S. Necco and R. Sommer, Nucl. Phys. B622 (2002) 328, hep-lat/0 108008. 
[17] M. Delia Morte and L. Giusti, in preparation. 



13 



